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ABSTRACT 



The random-phase-approximation (RPA) method with separableresidnal forces 
(SRPA) is proposed for the description of multipole electric oscillations of valence 
electrons in deformed alkali metal clusters. Both the deformed mean field and 
residual interaction are derived self-consistently from the Kohn-Sham functional. 
SRPA drastically simplifies the computational effort which is urgent if not deci- 
sive for deformed systems. The method is applied to the description of dipole, 
quadrupole and octupole plasmons in deformed sodium clusters of a moderate 
size. We demonstrate that, in clusters with the size > 50, Landau damping 
successfully competes with deformation splitting and even becomes decisive in 
forming the width and gross structure of the dipole plasmon. Besides, the plas- 
mon is generated by excitations from both ground state and shape isomers. In 
such clusters familiar experimental estimates for deformation splitting of dipole 
plasmon are useless. 



PACS: 31.15.Ew; 36.40.-c; 36.40.Cg; 36.40.Gk; 36.40.Vz; 36.40.Wa 



1 Introduction 



Most of free atomic clusters are known to exhibit quadrupole, hexadecapole, 
octupole, etc. deformations (see, e.g., and references therein). The defor- 



mations persist from hghtest |TT| to heavy |^ clusters and considerably influ- 
ence cluster's static and dynamical properties. For example, the axial (triaxial) 
quadrupole deformation leads to the splitting of the dipole plasmon into two 
(three) peaks p|, ^ and to specific low-energy orbital magnetic dipole reso- 
nance 0-111. 

Investigation of deformation effects in the dynamics of clusters is a topic of 
current interest, for reviews and monographs see |]T3|-|^. At the same time, the 
handling of excitations in deformed clusters is much more tedious because there 
are less selection rules and thus larger irreducible particle-hole spaces. In particu- 
lar, random-phase-approximation (RPA) calculations for multipole oscillations of 
the valence electrons (plasmons) require diagonalization of matrices with a very 
high rank (up to 10^ and more). 

For this reason, the calculations are mainly limited either to light clusters 
(with number of atoms N < 20) where the computational effort is still acceptable 
||TI| , or to very heavy clusters where one may neglect effects of quantum shells and 
use macroscopic models, thus avoiding a deal with a huge configuration space. 
The main computational troubles are concentrated in the intermediate size region, 
10 < < 10^, where sizable effects of quantum shells take place. Here there exits 
only few calculations with RPA or fully fledged Time-Dependent Local Density 
Approximation (TDLDA) (see, for example, |2^-[H). Thus one often chooses 
simplified methods, like Sum- Rule Approach (SRA) (see, e.g., |^) or local RPA 
(LRPA) p6|, They allow to reproduce the energy and gross structure of the 



plasmon, but cannot properly describe such important feature of cluster dynamics 
as Landau damping. At the same time, it is just Landau damping which greatly 
influences the shape of the plasmon and provides the dominant contribution to 
its width in medium and heavy clusters So, there is the need in RPA 

approach for deformed clusters, which, on the one hand, gives comprehensive and 
accurate description of multipole plasmons, including Landau damping effects, 
and, on the other hand, accomplishes this in economical way. 



To this end, a new kind of RPA method has been recently proposed [28, 3C 
It invokes a self- consistent expansion of the residual two-body interaction 
into a sum of separable terms: 



K 



{Pl,P2\Vres\hl,h2) l^kk'{Pl\Qk\hl){p2\Qk'\h2) ■ (l) 



k,k'=l 



Here, n^y is the matrix of strength constants, {p\Qk\h) is the particle-hole ma- 
trix element of the local one-body operator (5fc(r). The expansion index k labels 
various angular momentum channels A/i as well as a series of different radial func- 
tions within each channel. The method proposed is called as a "separable RPA" 
or SRPA. Due to the separable ansatz, SRPA has the advantage to transform the 



high-rank RPA matrix to a simple dispersion equation, which drastically reduces 
the computational effort. This is crucial for studies of the linear response in de- 
formed or large spherical systems, when one has to deal with a huge particle-hole 
configuration space. 

It should be emphasized that SRPA has important distinctions as compared 
with trivial separable RPA schemes widely used in many-body studies, especially 
in nuclear theory |^T], Such schemes exploit only one separable term with 
an intuitive guess for the separable one-body operator Q{r), and the strength 
constant k is fitted as to reproduce available experimental data. The present 
SRPA constitutes a systematic expansion of the self-consistent residual interac- 
tion. The analytical expressions for both strength matrix k^^/ and operators 
Qki'r) are uniquely prescribed by the given V^es- There are thus no adjustable pa- 
rameters. One can easily control the convergence of the expansion by comparing 
calculations with different cutoffs K. In fact, already a single SRPA separable 
term provides a fair description of the gross structures, and a good handful of 



terms allows to reproduce exact RPA calculations ||2^, The SRPA scheme 
was also proven to cooperate with the more involved Hamiltonians containing 
non-local pseudo-potentials 



The separable expansions p8[ and [BO] use different self-consistent prescrip- 



tions. The method deals with simple Qfc(i') operators but for the price of 



more separable terms. Vice versa, the operators in |^| are more involved but 



the method needs less number of separable terms. In the present paper, we will 
deal with the method [^] where the form of the separable operators (5fc(r) and 
strength matrix Kkk' are determined from the Vibrating Potential Model (VPM) 
[^-|^. The number of separable terms K is usually limited by 2-3 and 4-6 in 
light and heavy clusters, respectively. 

First SRPA results for the dipole plasmon and scissors mode in deformed 
clusters were presented in refs. |TB[-|T^, [^-[^. But these calculations were 



based on a simple shell-model ground state using a deformed Woods-Saxon po- 
tential where the deformations were not obtained within the same calculation 
scheme but taken from other sources, both experimental and theoretical. In the 
present paper, we present fully self-consistent SRPA scheme for deformed clus- 
ters. Both, the deformed mean field and residual interaction, are derived from 
the same Kohn-Sham functional. Ground state deformations are optimized by 
minimization of the total energy. SRPA is applied to detailed description of the 
dipole, quadrupole and octupole plasmons in axially-deformed singly-charged Na 
clusters of a different size. Both quadrupole and hexadecapole deformations are 
taken into account. 

We will show that in light clusters the dipole plasmon exhibits distinctive 
two-peak structure caused by the deformation splitting and just this splitting is 
mainly responsible for the plasmon width. The ratio between the peaks allows to 
conclude on either prolate or oblate shape takes place. The situation changes for 
medium and heavy deformed clusters where Landau damping becomes the main 
mechanism defining the width and gross structure of the plasmon. Moreover, 
such clusters exhibit a complicated energy surface with shape isomers which, in 



spite of a tiny energy deficit relative to the ground state, possess quite differ- 
ent deformation. Obviously, such isomers can contribute to the plasmon. Due 
to both the factors. Landau damping and contributions of shape isomers, gross 
structure of the dipole plasmon in medium and heavy deformed clusters cannot 
already serve as a direct fingerprint of the deformation and familiar experimen- 
tal prescriptions for getting magnitude and kind (prolate or oblate) of cluster's 
quadrupole deformation become useless. 

In Section 2, the deformed Kohn-Sham mean field is described. The existence 
of shape isomers in heavy clusters is demonstrated. SRPA equations are derived 
and specified for deformed clusters. In Sec. 3, SRPA results for the dipole 
plasmon are discussed. The same is done for quadrupole and octupole plasmons. 
The summary is given in Sec. 4. Appendixes A-C contain details of the formalism. 



2 Deformed Kohn-Sham mean field 
2.1 Basic constituents of description 

We start with the Kohn-Sham functional for a cluster including valence elec- 
trons and atoms 



where 



^tot(t) = ^kin(t) + E,,{t) + Ec{t) (2) 

E^Ut) = ^ j drT{r,t) , (3) 
Excit) = J drp{r,t)e^,{p{r,t)) , (4) 

Eait) = I /cirdr/^(^'^^-^'(y^^^^^^,^'^^-^-^^^^) (5) 
1 J J |r — Y\\ 

are the kinetic energy, exchange-correlation (in the local density approximation 
p8| , see Appendix A for details), and Coulomb terms, respectively. Here, Pt(r) is 
the ionic density in the jellium approximation, p(r, t) and r(r, t) are, respectively, 
the density and kinetic energy density of valence electrons, expressed through the 
single-particle wave functions ipqiv^t) as 

OCC occ 

p{r,t) = ^«;,|^,(r,t)|2 , r(r, t) = ^ u;,|V^,(r, t)|2 (6) 
<? g 

with a temperature T introduced through thermal occupation weights 

f^a = z — i; — ■ (7) 

' 1 + expC-^) 

Here, Eg is the energy of the single-particle state q. The chemical potential Xp 
serves to ensure the correct number of valence electrons in the static ground 



state, i.e. / drpQ^r) = N^. The temperature regulates the occupation numbers 
of the single-particle states, which becomes important for clusters with partially 
occupied subshells. In medium and heavy clusters where the energy surface is 
very complicated the introduction of a temperature softens numerical fluctuations 
and thus improves the convergence of numerical results. 

The Kohn-Sham single-particle Hamiltonian is obtained by variation 



This yields 



h{r,t) = -— v'+f/xc(r,t) + f/c(r,t), (9) 
2m 

UUr,t) = e.,(p(r,t))+p(r,t)^^:^^f^, (10) 

dp 

Ucir.t) = e / dri ■ ■ , (11) 



see Appendix A for more details. 

The ionic background is described in the soft jellium approximation. The 
density of the deformed jellium is 

^^^""^ " l + exp{{r-R{e))/a) ^^^^ 

where quadrupole and hexadecapole deformations, 62 and ^4, are introduced 
through cluster's radius as 

R{e) = i?o(i + E ^Ano(0)). (13) 

A=2,4 

Here Rq = CvsN^^^; Vg is the Wigner-Zeitz radius (r^ = 3.96 a.u. for Na clusters); 
the coefficient C is adjusted to ensure volume conservation / drpiiv) = N, piQ = 
is the bulk density. The diffuseness values a = 1 and 0.8 a.u. are chosen for 
calculations with high (T = 400 — 800 K) and low (T = 105 K) temperature, 
respectively. These values allow to describe the energy of the dipole plasmon in 



a wide size region ( see Ref. and present results). Decrease of the diffuseness 
parameter results in a redshift of the dipole plasmon which should follow lowering 
the temperature. The recent experimental data for the dipole plasmon in light 



and medium clusters ||39| were obtained at T = 105 K. So, just this temperature 
and the value a = 0.8 a.u. are used in the paper for all clusters except of the 
heaviest one Nafig. For the latter cluster the experimental data are absent and 
so the common value a = 1.0 a.u. is used. 



2.2 Representation of wave functions for stationary de- 
formed states 



Stationary states of the cluster are obtained as solution of the stationary Kohn- 
Sham equations hQ{r)ipg{r) = Sqipqir) where /io(r) is the single-particle Hamil- 
tonian in the static limit. We expand single-particle wave functions of the 
deformed Kohn-Sham potential in the complete basis of eigenfunctions 0s (r) = 
0niA(r) of the same but spherical Kohn-Sham potential: 

V'g=^n.A(r) = ^a^;0„iA(r). (14) 

nl 

The spherical states 0s (r) = 0nZA(r) = Rni{,r)Yip^(yL) are characterized by quan- 
tum numbers s = nlA (the node number, orbital moment and its projection to 
the symmetry axis, respectively), where the orbital moment is connected with 
the space parity as vr = (—1)'. For the sake of simplicity, spin functions are 
omitted though, of course, their contribution is taken into account in all relevant 
values. Deformed states are specified by asymptotic Nilsson-Clemenger quantum 
numbers q = Afn^A where A/" is the principle shell quantum number (the total 
number of quants, Af = rix + Uy + riz) and tt = (—1)'^. The levels in axially- 
deformed clusters are twofold (for A = 0) or fourfold (for A 7^ 0) degenerated. 
The electronic density reads in this representation 

occ 

Po(r) = 2 (2-5A,o)wg|V',(r)p. (15) 

Making use the multipole expansion for static potentials ([10|)-(|TT1) and densities 
([T^)-([T5|), the equations to determine single-particle energies Eg and wave func- 
tions ipq are derived. As a first step, the Kohn-Sham problem is solved in the 
spherical limit and then the basis obtained is used for deformed Kohn-Sham po- 
tential. At each iteration, the conservation of A'^e and A^ is controlled by fitting 
the parameters Xp and C. The procedure is performed for every point in the 
deformation map {62, 64}. The equilibrium deformation parameters 62 and 64 are 
determined by minimization of the total energy (^. Details of the calculation 
scheme are given in the Appendix A. 



2.3 Results for stationary states of deformed clusters 

We consider axially deformed clusters Naf^,, Na2Y, A^^ss, Na'^^, and Naiig. These 
clusters cover different size regions and represent both prolate and oblate shapes. 
Following calculations they do not exhibit neither triaxiality, nor octupole 

deformation and, so, we may safely deal only with axial quadrupole and hexade- 
capole deformations. The equilibrium values of 62 and 64 and the corresponding 
moments are presented in Table I. The moments are defined as 

An J drpoiryYxo 5 
u\ = = , K 

^ 3 NeR 



5 / drpo{r)r 
3 /rfrpo(r) 



A = 2,4 



As we checked, the calculated moments are close to those obtained within the 
Structure Averaged Jellium Model (SAJM) [0]. Table I shows that some clusters 
exhibit considerable hexadecapole deformation. Further, in particular clusters the 
values ^4 and /34 seem to be disentangled. For example, large ^4 corresponds to 
small Pi in Naf^ and vice versa in Na^j. This reflects the fact that hexadecapole 
moments are determined by both quadrupole and hexadecapole terms in (p^) 
and, moreover, the role of the quadrupole term is leading. 

The calculations show that in light clusters Naf^, Na2f and Na^^ the ground 
state deformations are characterized by a deep and distinctive minimum with 
considerable energy gain as compared with shape isomers. The number of isomers 
increases and their energy difference shrinks with increasing system size. The first 
isomers in Na'^^ and Naf ig are given in the Table I. They exhibit an opposite 
sign of quadrupole deformations as compared to the ground state. Fig.l shows, 
as an example, the energy surface of Naf ig, which has three distinct minima. The 
second isomer demonstrates a dominant hexadecapole deformation. 

In Figs. 2 and 3, the single-particle spectra of prolate A^a^y and oblate Na^^ 
are shown and compared with the corresponding spectra in the spherical limit. 
It is well visible how the deformation removes degeneracies. Note that prolate 
and oblate clusters exhibit the opposite assignment of quantum numbers inside 
subshells. 



3 Separable approach to dynamics 
3.1 Small amplitude dynamics 

Up to now, we have determined stationary states with a mean field hoi^r), den- 
sity po{r) and single-particle wave functions '?/'g(r). We are now proceeding the 
dynamical Kohn-Sham equations 

h{r,t)ij,{r,t) = t-^^^g{r,t) . (16) 

The many-body wave functions are composed from the single-particle wave func- 
tions ipq{r,t) as Slater determinants \E^(ri, fat, t). In the linear regime, when 
only small oscillations about the stationary state are considered, one gets in first 
order perturbation 

p{r,t) ~ po{r) + 6p{r,t), 
h{r,t) ~ ho{r) + Sh{r,t) 

where the response mean-field reads in detail 

6h{T,t) = (^)^=^„5p(r,t) + J = Tr^ {v;es,i25p2} • (17) 



The Sh{T, t) is a purely instantaneous and local operator by virtue of the local- 
density approximation inherent in the Kohn-Sham scheme used here. For har- 
monic oscillations it can be written as 



Further, for small-amplitude harmonic oscillations in the vicinity of the static 
ground state \l/o; the pertinent ansatz for a perturbation of the wave function 
reads 

5^>{t) = Y: «.e-* + c;,e— *) ala^^o (19) 



ph 

where aj, and ah are creation and annigilation operators of particle and hole states, 
respectively. Insertion of (|18|) and (p!9| ) into the time-dependent Kohn-Sham 
equation (p!6D with subsequent linearization and selection of the part oc e"""^* yields 
finally a standard RPA problem for frequencies uj and corresponding amplitudes 

{eph + <yuj)c;h+{p\5h''\h) = ^ (20) 

where a = ± and 

{p\5h^\h) = Y: {{ph'\V,,,\hp')c;'h' + {pp'\Vres\hh')c-,l,) . (21) 
p'h' 

3.2 Separable RPA 

The full RPA, although straightforward, can become very expensive because the 
dimension of the matrix of VJ-es grows huge in demanding applications. The 
separable approach (Q) is designed to reduce the expense of such calculations. 
Inserting (0) into the full response (pT]) yields 

{p\6h''\h) = Y^l{p\Q>^\h) , (22) 

k 

k' •' 

= Y.^kk'j:{iP\Qk'\hKh+{h\Qk'\pKH^). (23) 

k' ph 

Note that the values 6h'^{T^, a% and 5p"{v) become independent of a for purely 
local separable operators Qkij^)- Such situation takes place when the motion is 
driven only by variations of time-even densities. This is just the case for clusters 
where the dynamics is determined by oscillation of time-even density of valence 
electrons. Then, 

ak = Y. i^kk' Y.^p\Qk'\h) (c+^ + c-^ . (24) 

k' ph 

Using (|22|) , the RPA equations (pOD can be reshuffled into 



± _ ^ llk<^k{p\Qk\h) , . 



Inserting P3| ) into (^) yields finally a system of linear homogenous equation for 
weights ak- 

K 

Skk'{uj)ak' = (26) 

fc'=i 

where 

f ^ ^<p\Qk\h><p\Qk'\h> e.ph 1 . . 

Sfcfc'(a;) = 2^ -2 — — . (27) 

In Eq. (|27|) the sum runs over all particle-hole states of a given multipolarity. 
The condition 

det I s{u) 1= (28) 

gives SRPA dispersion equation for eigenvalues uj. 

The rank of the SRPA matrix (^) is equal to the number K of the separable 
operators, which, as is shown below, is typically 3-6. So the rank is dramatically 
less than in the case of involved RPA methods. This results in drastic simplifi- 
cation of RPA calculations. At the same time, the total number of SRPA roots 
equals to the number of particle-hole configurations and every SRPA state is a 
superposition of these configurations. The SRPA keeps the important RPA fea- 
ture to describe the Landau damping but gets this with much less computational 
effort. 

It's straightforward to show that Eqs. (^6])-(^) can be also obtained with a 
generalized schematic RPA Hamiltonian 

H = Ho-^J2 ~^kk'QiQk' ■ (29) 

^ kk' 

This can be done by substituting the Hamiltonian (|29D to the standard RPA 
equations 

[H,d]]=u;,d], [H,Q,] = -LoA, [0,,0j,]=5,y (30) 



for the creation phonon operator O] = Y^ph y^ph^p^h + ^ph^h%) hermitian 
conjugate. The matrix [kkk'] in (|29|) is inverse to the matrix [n^k'] (H 



''kk' 

Note that relative contributions of operators Qk{T^) to RPA states are not 
fitted but regulated self-consistently by the weights dl- Every RPA state j has 
its own set of the weights. The normalization condition for the particle-hole 
amplitudes c^^ and weights al can be written through derivatives of Skk'i^j)'- 



ph kk' ph \^ph ^j) 



1 _ ■ _ ■ d 
ll^Wk'-^Skk'{.ujj) = 2. (31) 



3.3 Self-consistent prescription for separable operators 

Up to now, the operators Qfc(r) and strength matrix Hkk' were not still chosen. 
Obviously, the success of SRPA depends on the choice. The better it is, the 
smaller number of the separable operators we need to reproduce V^es- We are 
considering systems with pronounced collective modes (the various plasmons in 
case of clusters) and so just collective deformations should motivate the form of 



Qk{^)- Here we will follow the VPM prescription |]3J] which was very effective in 



the SRPA calculations for spherical clusters [^. In the VPM, the operators Qk{^) 
are self-consistently generated by collective variations of the density which in turn 
is determined by the scaling transformation P^, HT|. Being thus constructed. 



only few operators are enough to cover successfully the dynamics. Specifically, 
the collective deformation is achieved by the scaling transformation 

K ^ 

^(t) = H e"fcWEJi['^o{r,),A.(r,)]^^ _ (32) 

k=l 

Here, the local operators /^(r) determine the intended time-even shape of the 
deformation. For example, f = z leads to a translation along ^-direction, / = 
r^y2o results in an ellipsoidal deformation. The collective amplitudes ak{t) = 
akcos ujt determine the relative weights of different deformation modes. For small 
deformations otk-, the density variation becomes 



5p(r,t) = {^>{t)\p\^>{t))-{^o\m 



oj 



K 



- E«feW(V/Oo(r)- v/fe(r) + Po(r)A/,(r)) (33) 

k=l 

which, after submission into (|l^, results in the response Hamiltonian 

K 

Sh{r) = J2akQk{r) (34) 



k=l 

composed from the operators 

Qfe(r) = J driKos(r - ri) (vPo(ri) • V/fc(ri) + po(ri) A/fe(ri)) (35) 

with 



2 



Kes(r - rO = (-^)p=po<5(r - n) + . (36) 

op |r — ri| 

Operators (|35| ) take, by construction, into account all deformation distortions of 
the ground state (this point is discussed in the Appendix B where the explicit 
expression for Qkij) in deformed clusters is given.) It is now self-suggesting that 
these operators are identical with the basis operators to be used in the separable 
ansatz This amounts to identify the collective response ( plj) of operators 
([35| ) with that from the separable ansatz (p2D, which finally yields the strengths 
coefficients as 

K,l, = -J drQ,(r)(vPo(r) ■ v/fc'W + Po(r) A/,,(r)) . (37) 



Eqs. and (P3[)-(P7D constitute the complete set of self- consistent 

SRPA equations. One has to point out that, hke involved RPA, SRPA maintains 
full details of Iph structure in the spectra. The only point is that the residual 
interaction is parametrized in terms of collective picture. 

It is easy to see from (|3^) that just input local operators fk(j) determine the 
character of the perturbation. The choice of the operators is the most delicate 
point of the approach and devotes a special attention (see discussion in Appendix 
C) . We use the set of hermitian operators: 

/A.P..(r) = r'''^''{Y,^,{n) + Yl^m (38) 

with 

XkPk = 10, 12, 14, 30, 32 for (/x = 0, 1), 
XkPk = 20, 22, 24, 40, 42 for (/i = 0, 1, 2), (39) 
XkPk = 30, 32, 34, 50, 52 for (/x = 0, 1, 2, 3) 

for dipole, quadrupole and octupole excitations, respectively. It is seen that in 
all the sets the first operator has the form of the applied electric external field 
in the long-wave approximation, the others two represent a variety of the radial 
dependencies to cover different slices of the system and the last ones take into 
account the coupling of different A-modes with equal /i and parity, which takes 
place in deformed systems. Relative contributions of different input operators 
fxkPkiJ.{^) to residual interaction are self-consistently regulated by the strict 
computation of the strengths according to Eq. (^). 



3.4 Response to electric fields 

The photoabsorption cross section for excitation of the state j = in axially 
deformed cluster has the form 

a{E\^^■ gr ^ j) = ^^^^^ + V^a.. (40) 

where 

Ma^, = <j\r^Yx^\0> (41) 

= 4^ E < P\r'Yx,\h > (4+ + c^p-) = Yl (^IM^j) 

is the reduced matrix element of EXfi transition from the ground state to the 
RPA state j and 

. / X Eph < p\r^Yxf,\h >< p\Qk\h > 

^k[UJj) = 2^ -2 -2 . (42j 

ph ^Ph - 

Plasmons in deformed clusters involve a lot of RPA states which in any case 
cannot be well resolved. So, it worth to consider not every RPA state but an aver- 
aged response. Strength function formalism allows to get the averaging response 



without dealing with RPA problem for every state, which gives a huge gain in the 
computational effort. The factorized residual interaction used in SRPA makes 
derivation of the strength function especially simple |^ . Following this line, the 
strength function is defined as 

criEXfi, iu) = Yl cr(£'A/i; gr t-^ j)T]{LU - Uj) (43) 
j 

where 

^^"-"^•^ = ^ (..-c./+(A/2)^ ^''^ 

is Lorentz weight with an averaging parameter A. Using Cauchy residue theorem, 
one may rewrite (^31) in the form independent of j (when the response is directly 
computed for given u) . The derivation starts with the squared matrix element of 
EXu transition which, provided that SRPA matrix Skk' (Eq- (pT])) is symmetric, 
can be written as 



- ^^^^ 

where 



= J2^kk'{uJj)Ak{ujj)Ak'{ujj), (46) 

kk' 

and Skk'i'^j) is the algebraic supplement of the matrix element Skk'i'^j)- Then, 
provided that the function 1/ det |s| has only one-multiple poles u = ±ujj, the 
strength function ( |i3| ) is expressed through the corresponding residues on the 
complex plain z: 

a{EXf,,u) = 8^ %^(2A 1 1)!!]^ (M^^-^ E ^^^[-^(^^ ^)]^=±^. (47) 



where 



det \s{z)\ 



Following Cauchy theorem, sum of all the residues (covering all possible poles 
of the function (^§1) ) is zero and so we can express the residues with z = Uj 
through the rest of the others: 

Res[Fl=^^= - {Res[Fl=_^^ + Res[Fl^^ 

+ i?es[F],=,±i(A/2) +i?es[F],=±,^J . (49) 

It's easy to prove that lim|^|_^oo -^(^? -z) = for A = 1,2. Also, Res[F]z=-u;j and 
Res[F]z=-£pf^ can be neglected for large positive z- values in the plasmon energy 
region and for relevant values of the averaging parameter A. Remaining residues 
over the poles z = u ± i(A/2) and z = Sph give the final outcome 

a(£A/,,.) = 8,»^p^±l^^^ (50) 



X 



-$5 



TT 



Z 



2A-1 



G{z) 



det 



=a;+j(A/2) p/i 



Unlike expression for the standard strength function ||3T 
from two terms 



ig, Eq. (|0D consists 



The first one is contribution from the residual interaction and 
the second one is unperturbed strength function. This form is convenient for the 
analysis. The standard strength function can be reduced to the form (|50|) if the 
width A is introduced to ± ti;j)-distribution. This gives the Lorentz shape. 
Eq. (|50D can be used for calculation of El and E2 photoabsorption. 
Model-independent energy-weighted sum rules (EWSR) 



1^ ph Mi 



2 



-A(2A + lyNe < r'^'' > 



(51) 



is used to control the completeness of the Iph configuration space. 



4 Results and discussion 
4.1 Dipole plasmon 

Results of SRPA calculations for the dipole plasmon are exhibited in Figs. 4-6. 
Two different presentations are used. First, to demonstrate the complex structure 
of the plasmon and Landau damping, the photoabsorption strength (^) is given 
for every RPA state by a vertical bar (in eVA). Second, for the convenience of 
the comparison with experimental data, the strength function (^) averaged with 
the Lorentz weight is used. The value of the averaging parameter is A = 0.25 eV 
to simulate the typical thermal broadening of the plasmon. 

In Fig. 4, the photoabsorption is presented for three light deformed clusters, 
prolate Naf^ and Na2j and oblate Na^^^. The SRPA calculations are compared 
with recent experimental data |^ obtained at rather low temperature T = 105K. 



The agreement with experiment is very satisfying. It should be mentioned that 
the average energy of the plasmon was slightly corrected by decrease of the dif- 
fuseness parameter a in (|T^) from 1.0 to 0.8 ao. The value 1.0 ao chosen to 
simulate ionic effects at high temperatures (T = 400 — 800K) |j2^ is not suitable 
for T = 105-?^. The decrease of the surface diffuseness a leads to a blueshift of 
the plasmon, which is indeed seen when lowering the temperature. Beside the 
good description of the average energy, the deformation splitting of the plasmon 
into two peaks is well reproduced. The peaks correspond to /i = and 1 dipole 
modes. The latter involves contributions from both = ±1 projections and so is 
about twice as large. This mode forms the upper part of the plasmon spectrum 
in prolate clusters and the lower part in oblate ones (see also Fig. 9 e)). The 
splitting is well developed in the prolate clusters and less strong in the oblate 
iVogg. Though first evidence of the Landau damping is already seen, the plas- 



mon gross-structure and width are still mainly determined by the deformation 
splitting. 

Following such a splitting analysis, one is tempted to associate with A^a53 (see 
Fig. 5) an oblate deformation. Its plasmon shape is rather similar to the one for 
oblate Na^c^ and the broad peak with a small right shoulder at 3.1 eV can, in 
principle, be approximated by two Lorentz curves "justifying" an oblate shape, 
as done e.g. in However, our calculations show that this interpretation is 



misleading. Unlike lighter clusters, Naf^ shows considerable Landau damping 
(see bars in Fig. 5, top panel). This broadens the plasmon and modifies its gross- 
structure. In particular, the structure at 3.1 eV is caused not by /i = mode 
but by the group of yU = 1-states. Our calculations predict for this cluster the 
prolate shape, which results in very satisfactory agreement with the experimental 
data. Good description is also obtained for other clusters in this size region 
(see the comprehensive SRPA analysis of the data in |^3|). Our conclusion 



on prolate shape of clusters in this size region agrees with calculations of other 



groups Altogether this means, that generally accepted experimental practice 
to treat plasmon shape within simple deformation models and thus to conclude on 
cluster deformation is not appropriate for medium-size clusters. This way does 
not take into account important role of Landau damping in forming plasmon 
gross-structure and width and so can be misleading. 

For medium and large clusters, the disentangling of the dipole spectrum be- 
comes even more complicated because of possible incoherent contributions from 
isomers. Table I shows that A^a^3 has oblate isomeric state with a tiny energy 
deficit as compared with the ground state. This isomer can also contribute to the 
observed spectrum. The bottom panel of Fig. 5 demonstrates that the isomeric 
state yields similar pattern, in spite of reverse ordering of /i = and 1 modes. In 
A^a^3 both the modes are broad and have two-peak structure. In some important 
details the isomer plasmon deviates from the experimental data. 

The picture is complicated with further growing cluster size. As was shown in 
Fig.l, the energy surface in Nafig has already three distinctive minima: oblate 
ground state, first prolate isomer and second hexadecapole isomer. The corre- 
sponding plasmons are exhibited in Fig. 6. Contributions of /i = and 1 modes 
are also given to show the role of the deformation. All the plasmons demonstrate 
strong Landau damping which determines, in a large extent, their width. Just be- 
cause of Landau damping, the plasmons, in spite of different cluster deformations, 
show very similar broad shapes. Obviously, all three deformation configurations 
of A^ctiig can contribute to the measured spectrum. Indeed, clusters exhibit ther- 
mal shape fluctuations which are incoherently added in experiment. The shapes 
of the ground state and first isomers have to provide the dominant contribu- 
tion. The observed plasmon reflects a statistical mix of many shapes with a few 
dominants. 

To summarize, the calculations demonstrate increasing importance of Landau 
damping and isomer mixing for the dipole spectrum. Already at medium cluster 
sizes these factors cannot be neglected. As a result, the often used prescription 
to conclude on cluster shape on the grounds of approximation of plasmon shape 



by 1-3 Lorentz curves is insufficient for these system sizes. 

The strong fragmentation of the spectra through Landau damping may raise 
doubts on the collectivity of the plasmon oscillation. The strength is distributed 
over many sub-states and calculations show that even the strongest RPA states 
do not exceed more than 2-4 times the strength of particle-hole configurations, in 
heavy clusters even less. Nevertheless, the entity of RPA states near the plasmon 
frequency can still be called collective. One may view this as a two step process. 
The dipole excitation ffist generates a collective surface plasmon state which 
is then quickly distributed over the nearby particle-hole states through Landau 
damping. The lifetime of the collective dipole excitation can be estimated to be 



about 10 fs in large Na clusters |21]. The collectivity is also testified by reach 



particle-hole structure of RPA states. 

4.2 Quadrupole and octupole plasmons 

In Figs. 7-8 the photoabsorption cross sections for quadrupole and octupole 
plasmons in prolate Na2-j and oblate Na^^^ are presented. Contributions from all 
the projections /i are given to demonstrate the role of the deformation broadening. 

The figures show strong Landau damping. Its contribution to the width and 
gross structure of plasmons is of the same scale or even larger than from deforma- 
tion splitting. The prolate Na^-j has a strong deformation and we see a distinctive 
blueshift of /i-branches of the plasmon with increasing /i. This makes the defor- 
mation splitting comparable with Landau damping. The deformation is weaker 
in the oblate Na^^, which makes the deformation splitting less pronounced than 
Landau damping. 

Both, quadrupole and octupole plasmons, are well concentrated (especially 
in the less deformed Nat,^), which helps their experimental observation. The 
plasmons have a high-energy tail. Analysis of structure of main RPA states 
from the plasmon region shows that they have about the same character as in the 
dipole case. There is initially a collective multipole plasmon which is then quickly 
Landau damped and fragmented over the many nearby particle-hole states. 

As is seen from the figures, photoabsorption cross sections for quadrupole and 
octupole plasmons are 10^ and 10^° times weaker than in the dipole case. Both 
plasmons share the energy region with the all dominating dipole plasmon and 
so should be completely masked by the latter. Strong Landau damping of the 
plasmons additionally decreases their chance to be observed in photoabsorption. 



However, following the estimates inelastic electron scattering could resolve 
these states. At certain angles of outcoming electron, contributions of A = 2 and 
3 plasmons to the cross section can successfully compete with the dipole contri- 
bution. The reaction is very demanding as it requires an intense, monochromatic 
and well coUimated low-energy electron beam. However, progress of experimental 
techniques gives us some hope. SRPA results can be used in such experiments as 
a ffist guide for plasmon energies and strength distributions. 



5 Conclusions 



The separable RPA (SRPA), a self-consistent approach to RPA, is presented in 
connection with investigation of multipolc plasmons in deformed clusters. Both 
ground state properties and RPA response are calculated using the same Kohn- 
Sham functional for axial shapes. Due to self-consistent factorization of the resid- 
ual interaction, SRPA dramatically simplifies RPA calculations. At the same 
time, the method ensures high accuracy of the calculations. As compared with 
other RPA methods, SRPA provides a comprehensive treatment of Landau damp- 
ing with a minimal computational effort. This is crucial for deformed systems 
where one deals with an impressive particle-hole configuration space. 

SRPA has been applied to dipolc, quadrupole and octupole plasmons in a va- 
riety of axial clusters. The main attention was paid to the dipole plasmon. Light, 
medium and heavy clusters of both prolate and oblate shapes were covered. Both 
quadrupole and hexadecapole deformations were taken into account. Contribu- 
tions to the plasmon from both ground and isomeric states were analyzed. Good 
agreement with available experimental data was achieved. 

Unlike light axial clusters, where dipole plasmon has a typical two-peak struc- 
ture reflecting the deformation splitting, the dipole spectra in medium and heavy 
deformed clusters show usually a broad bump with relatively small variations of 
the shape. The deformation splitting is masked by strong Landau damping which 
inhibits a direct determination of ground state deformations from the spectra. 
Moreover, such clusters have shape isomers which further overlay the spectra. As 
a result, a correct treatment of plasmon properties is possible only on the grounds 
of fully microscopic calculations taking into account all these effects. SRPA is a 
promising method for STich investigations. 

Quadrupole and octupole spectra were shown to stay well concentrated in 
spite of the deformation splitting and Landau damping. Both spectra have a 
considerable high-energy tail. The quadrupole and octupole spectra cannot be 
observed directly in photoabsorption experiments but there is some hope to dis- 
entangle them in inelastic electron scattering. 



Appendix A 

The problem (y) for the deformed Kohn-Sham potential 



2m 



V' + UUr))ijg = Eg^Pg (52) 



can be rewritten as a system of equations for amplitudes of the expansion 
(0): 

{E, - eg)ai + < <Ps\ifcir) + U.c{r)\(ps' >= (53) 

<?' 

where 

f/c(r) = f/c(r)-f/0(r), (54) 
Mr) = f/.c(r)-f/°,(r), (55) 
and values U^{r) and U^ci'"^) enter the problem for spherical Kohn-Sham potential 

{-^ +U'c{r) + Ul{v))<P, = E4s. (56) 



For the quantities ([T^), (p!5|), (^4]), and (|55|) the multipole expansion is used. In 
particular, the electron ground state density ( plSf ) is written as (see notations in 
Sec. 2) 

Po(r) = E>3vo(^)PA(0, (57) 



with 



paW = I . E (-1)^(2 -M^. 

y 7r(2A + 1) q=Mn^A 



X E V(2^2 + l)(2/i + l)Qto,u,o<-A,/,A (5J 



"•2,^2, "1,^1 



X a. 



n2l2 

In all the multipole expansions, including Eq. (|57D we take A = 0, 2, 4, 6, 8, which 
provides sufficient numerical accuracy. 

Kohn-Sham equations (p3D together with the density equation ( |57D are solved 
by an iteration method. Eigenvalues Eg and eigenfunctions (ps for the spherical 
Kohn-Sham potential are used as the input. 

The equilibrium deformations S2 and ^4 are determined by minimization of 
the total energy (^. It is convenient to use for the kinetic energy the expression 

1 2 occ 

Et = — (2-MEE«^ (59) 

X ( J drr^ — — — — + /(/ + 1) y drRn^i[r)Rn^i[r)). 



where (-1)' = tt, A > 0. Further pS 



with 
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Ry ao, r,(po) 
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-0.066i??/{(l + x^)log{l + x-^) + -x 



11.4 ao 

The Coulomb term can be written as 

o2 



Ec = ^J dr{po{r) - p,{r))Uc{r) 



Appendix B 

Here we present explicit expression for the operator of residual interaction 
(|35| ) in deformed clusters. Making use a multipole expansions for densities (|l^) 
and (|TB|) and potentials (Q) and (|5BD, one gets 

Qk{r) = E T.i'^L.m + yLm m 

A L 

x[^A'ol.Qi?A(-) + EE(>^-.(^) + >^LV(fi))<'l,QiltAv(0] 
A' 



where 



47re^ 



^'''^^ y47r(2L + l)2L + l 

r 00 

X |r-''+'' / dnn'l^l{ri)ri" + r'- j (;riK|.'^(ri)rr''""] 



Here, Px('") ^^^{r) are expansion coefficients in (|58|) and U^cij) = Y^x f^i'c(^)^Ao(^)' 
respectively, and 



Kl = (p. + A. + 1)A(2A.-1)(5S.-4a.^ 

- (Afc - Pfe) V(Afc + 1)(2A, + 3 )(i?g, - Agj , 

M^S = (P^ + A. + 1)v/A,(2A, - 1)(( A + + 

+ (A, - p,)^{X, + 1)(2A, + 3)((A + + Mgj 



+ (Afc - R.) (pfc + Afc + l)v/(2Afc + l)(2A + l)C; 



LO 

AOAfeO 



where 



4L = V(AH-l)(2AH-3)(\+/ ^^^-1 t)c'".o^.-.o- 



LO 

A+1 Xk+1 ' 



LO 

A-1 Afc+1 • 



The values Cfi; , and those in big parentheses are Clebsch-Gordan coefficients 
and 6j-symbols, respectively. 

Expression ( |63D shows that coupling of oscillations of given multipolarity (Afcyu) 
with the spherical (A = 0) and deformed (A = 2, 4, ..) parts of the density leads to 
a family of terms with \ X — Xk \< L < A + A^. Thus the residual interaction driven 
by the operator (^) takes into account not only the deformation distortions of 
the ground state density but also the coupling of /i'^ modes with different A. This 
coupling is pertinent only to deformed systems (see, e.g. discussion in [pTSf). 



Appendix C 

Here we will comment some essential points concerning the scaling (|3^), choice 
of input operators /^^.^^.^(r), and convergence of the SRPA results. 

Time-odd hermitian generators Px^^pi^^^r) = —i[hQ{r), fx^p,,^{r)] entering the 
scaling result in time-even paths ak{t). The generators can be considered as 
momenta conjugate to time-even hermitian coordinate-like operators fx^p^^iir). 
Since in our case all the dynamics is provided by the time-even electron density, 
there is no any need in time-even generators. This explains why the scaling 
includes not full but single-particle Hamiltonian. Indeed, if separable residual 
interaction in a schematic full Hamiltonian involves only time-even hermitian 
Qfc(i')-operators, then 

[HJ]ph^[hoJ] - Y.^kk'{Qk<^o\[Qk'J]\^o> 

kk' 

+ <^o\[Qk,mo>Qk') = [hoJ], (64) 

since averaged commutators between hermitian operators of the same time-parity 
is zero. 

It is easy to see from ( P^j ) that collective motion is determined by input 
operators f\t.p^^{r). The choice of the operators is motivated by physical reasons. 
The number of the operators should be as small as possible to minimize the 
computational effort and, at the same time, should be sufficient to ensure the 
convergence of the results. Specifically, if we study response of the system to an 
external field, then one of the operators must take the form of such field. One 
input operator is usually not enough to reproduce the true residual interaction, 
even if the separable interaction is constructed self-consistently. To overcome this 



trouble, we exploit the idea of the local RPA p6|, ^ to use not one but several 
input operators. In the local RPA the cluster is treated as a system of several 
coupled oscillators, which allows to describe the gross structure of the resonance. 
Instead, SRPA can be viewed as a system of coupled schematic RPAs. The 
operators fxf.pf.^ii'r) are chosen so that the corresponding operators of the residual 
interaction, Qk, have maxima at different slices of the system, including both 
surface and interior. Such input operators can be intuitively associated with 
different external probes causing the response of different parts of the system. 
Besides, several operators allow the system to relax additionally, which finally 
leads to high accuracy of the method ||29|. The number of input operators is 



dictated by convergence of the results and usually 3-6 operators are enough. 

Fig. 9 demonstrates the relative role of different input operators in forming the 
dipole plasmon. Prolate Na27 oblate Na'^^ are considered as typical examples. 
It is seen that the residual interaction posed by the leading input operator rYi^^ 
shifts the unperturbed dipole strength from 1-1.5 eV to the energy ~ 2.5 eV. 
This operator is enough to get the correct excitation energy. However, there still 
exists an artificial high energy strength at 3.5-5 eV, which is confirmed neither 
by experimental data nor by the TDLDA results. In fact it is a consequence of a 
poor separable approximation. Involving additional dipole input operators (panel 



c)) weakens the artificial strength. At the same time, the first operator is stiU 
leading. Its contribution remains to be overwhelming in most of the RPA states. 
The panel d) demonstrates the minor role of octupole operators responsible for 
the deformation coupling of the dipole and octupole modes with /j, = 0, 1. Our 
calculations performed for dipole, quadrupole and octupole plasmons in a variety 
of deformed clusters show that the deformation coupling of different multipole 
modes with a given fi is usually negligible. However, the deformation leads to a 
strong splitting of the plasmons to branches with different fi (see panel e) and 
Figs. 3-6). The comparison with experimental data for the dipole plasmon, given 
in Figs. 4-5, shows that SRPA correctly describes the deformation splitting. 

Fig. 10 compares the dipole strength calculated within exact RPA and SRPA. 
Cluster A^a53, being in focus of the discussion above, is chosen as a typical example 
for the comparison. For the simplicity, only /i = dipole mode is considered. The 
figure demonstrates good agreement which justifies high accuracy of SRPA. 
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Figure captures 



Figure 1: Energy surface (in Ry) of Nafig as a function of quadrupole and 
liexadecapole deformations. Tlie iso-energy counturs are separated by intervals 
0.01 Ry. Three distinct minima are marked by crosses (the energies at the minima 
are given in Table I). 

Figure 2: Single-particle spectrum in cluster Na2i ^^^o equilibrium (pro- 
late) deformations. In the latter case the Fermi level is A/'n^A = 321. The 
evolution of the spectrum with deformation is marked by dashed lines. 

Figure 3: The same as in Fig. 2 for oblate cluster Na^j- The Fermi level is 
Mn^k = 310. 

Figure 4: Photoabsorption cross section for dipole plasmon in small deformed 
sodium clusters. The deformation parameters are given in boxes. The experimen- 
tal data (triangles) are compared with SRPA results given as bars for RPA 
states (1^) and as the strength function ( ^3]) smoothed by the Lorentz weight. 
The bars are given in eVa.o. 

Figure 5: Photoabsorption cross section for dipole plasmon in Na'^^. The ex- 
perimental data (triangles) are compared with SRPA results for excitations 



from the prolate ground state (top) and oblate isomer (bottom). Contributions 
to the strength function from /i = and 1 projections (the latter has twice larger 
strength) are given by dashed curves. See Fig. 4 for notations. 

Figure 6: Calculated photoabsorption cross section for dipole plasmon in Nafig. 
The excitations from the oblate ground state (top), first prolate isomer (middle) 
and second hexadecapole isomer (bottom) are considered. See Figs. 4 and 5 for 
notations. 

Figure 7: Calculated photoabsorption cross section for quadrupole plasmon in 
prolate Na27 oblate Naf^. Contributions of projections with /i = 0, 1, 2 and 
their sum (bottom panel) are given as bars (|40|) for RPA states and as strength 
function (|43|) . The bars are given in eVaQ. 

Figure 8: The same as in Fig. 7 for octupole plasmon. 

Figure 9: Calculated photoabsorption cross section for dipole plasmon in prolate 
Na2j and oblate Na^^ (see notations in Fig. 4.). The cases include: a) unper- 
turbed (without the residual interaction) plasmon; b) with residual interaction 
generating by leading input operator Aipi = 10 : fiQ^i{r) = r(Yi^(f2) + Y'^j(f2)); 
c) with operators X^Pk = 10, 12, 14; d) with operators XkPk = 10, 12, 14, 30, 32, 
where two last ones provide the coupling with octupole modes; e) contributions 
of /i = (dashed curve) and /i = 1 (solid curve) projections to the cross section. 

Figure 10: Comparison of the dipole strength with /i = (as a squared transi- 
tion matrix element (PD) calculated within exact RPA (solid curve) and SRPA 
(dashed curve) in prolate Naf^. 



Table 1: Deformation parameters 82 and ^4 (as defined in Eq. (|TB|)) and 
quadrupole and hexadecapole moments [32 and [3^ (|16D , calculated for the ground 
and isomeric (in A^a^a and A^a^j^g) states. For isomers the energy deficits AE are 
given. 
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